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1 Schematic energy level for a chromophore. The electronic states are repre¬ 
sented by solid horizontal lines and the vibronic states are represented by 
dotted horizontal lines. S{ represents a singlet state and T* represents a 
triplet state. The photon absorption excitations are represented by solid 
vertical lines, and the decay processes are represented by wavy lines. The 
absorption coefficients and the decay constants kij are described in the 
text. The physical values used for our calculations are: < 7 0 i = 2.4 x 10 - 18 cm 2 , 
0 \2 = 3.0 x 10 - 17 cm 2 , (734 = 4.8 x 10 - 17 cm 2 , k\o = 0.144ns -1 , & 2 i = l.Ops -1 , 
k \3 = 77.Sms -1 , &30 = 50.0ms -1 , k i3 = l.Ops -1 , and k n , ^ 22 , k 33 are due to 


vibrational decays which are assumed to be instantaneous. 19 

2 Space-Time grid for pulse propagation. Squares represent coordinates where 

pulse is calculated and circles represent coordinates where carrier densities are 
calculated. 20 

3 Analytic and numerical solution for the fundamental soliton propagation. The 


solid line shows the analytic solution and the circles show the numerical solution. 20 
4 Relative error of and numerical solution for the fundamental soliton propagation. 21 
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Numerical study of short optical pulse 
propagation in nonlinear reverse 
saturable absorbers 

1 INTRODUCTION 

Recent developments in short pulse laser systems have introduced growing studies in the 
problem of ultra short pulse propagation in nonlinear media. Competitions among the effects 
of nonlinearity, diffraction, and group velocity dispersion (GVD) give rich spatiotemporal 
phenomena such as self-focusing[l], self-defocusing, and pulse splitting. Understanding of 
short pulse propagation is very important for protection of human eyes and optical equip¬ 
ment from high power laser pulses. Optical limiters are needed that give high transmittance 
at low input intensities and low transmittance at high input intensities. There are vari¬ 
ous mechanisms that result from nonlinear optical responses such as nonlinear absorption, 
nonlinear scattering and nonlinear refraction[2], Desirable properties for optical limiting 
candidates include, a high linear transmittance, a potentially low limiting threshold, a rapid 
response of picoseconds or faster, a broad spectral response, and a large dynamic range. 
Some semiconductors, gases and liquid crystals are used as optical limiting materials[3, 4], 
Liquid crystals exhibit high refractive nonlinearities but low response time (nanoseconds). 
Chromophores exhibiting nonlinear absorption, such as reverse saturable absorption (RSA), 
are under consideration for optical limiting applications^, 6, 7]. In a molecular system, 
RSA arises when the excited state absorption cross section is larger than the ground state 
absorbtion cross section. The process is modeled by several vibronically broadened elec¬ 
tronic energy levels. The general situation involves a five level system. The energy states 
included are three levels of the singlet state coupled to two levels of an excited triplet state. 
The mechanism of RSA is described in terms of simple rate equations coupled to a one¬ 
dimensional propagation equation of the optical pulse intensity. Figure 1 shows the energy 
level diagram for the five-state model. Other materials and mechanisms have been used such 
as fullerenes[8, 9], carbon black suspension[10,11], and two-photon absorption[12, 13]. While 
the main material studied in this paper is RSA, these materials are usually incorporated in 
a host material. With the recent introduction of novel high power lasers, the propagation 
of light in nonlinear materials is often subject to various phenomena such as self-focusing or 
self-phase modulation[14, 15, 16, 17, 18]. These effects are just beginning to be studied in 
RSA materials. Therefore, in order to study these effects we derived an equation[19] for the 
propagation of the electromagnetic field coupled to the rate equations because the traditional 
equations, which use only the pulse intensity, cannot answer these more detailed questions. 
Yet this equation cannot be solved analytically, therefore in this report we use numerical 
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methods. The principle numerical method is the split-step beam propagation method (BPM) 
that is coupled to a finite difference technique. One significant result is that, under certain 
circumstances, self-focusing can enhance the power limiting properties of RSA materials. 


2 NONLINEAR WAVE EQUATION 

From Maxwell’s equations we obtain the following wave equation 

( v “ &di?) E ^ r ^~ = 0, ^ 

where we assume that V • E = 0 and E || P. The electric polarization can be expressed as 
a sum of two parts such as the linear part Pi and the nonlinear part Pni 

P(r, t) = Pi{r ; t) + P n i (r, <), (2) 

and the linear part is related to the electric field by the following relation 

Pi(r, t) = e f X (t - t') ■ E(f, t')dt'. (3) 

Then the wave equation in the frequency domain becomes 

v2 + e ( w )“T E{r,u) + ^- 1 P nl (f,u) = 0, (4) 

C J CqC 

where 


e(w) = 1 + xM, 

_ TOO 

£(f» = / E(r, t)e~ iuJt dt, 

J-oc 

fOO 

Pm(r,co) = / Pmif^e-^dt. 

J -00 


In the slowly varying envelope approximation, it is useful to separate the rapidly varying 
part of the electric field as 

E(r,t) = A(r,t)e ikoz - iu,ot + c.c., (8) 

P nl (f, t ) = B(r, + c.c, (9) 

and A and B are slowly varying envelopes which satisfies the following conditions 

« *oMI. ^ (10) 

<oj„\A\, (T <<ul0 \Bl (11) 
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where Xi x, y , z. In the frequency domain, the electric field and the nonlinear polarization 
become 


E(r,u) = A(r,u-u) 0 )e +lkoz + A*(r,oj + Lu 0 )e~ ikoZ , (12) 

Pni{r,uj) = B(r,u-u 0 )e +ikoz + B*(r,w + u 0 )e~ ik ° z . (13) 

Then the wave equation becomes 

uj 2 

A(r,u-w 0 ) + — -B(r,u-u) 0 ) = 0, (14) 

6qC 


v 2 + 


' d V . m 2 

Yz +,h ) +e(a,) ^ 


where = d 2 /dx 2 + d 2 /dy 2 and we assume that e l{kz ^ terms and e~ i(kz ~^ terms are 
separable. Defining k(u) — yje{uj)uj/c we get the Taylor series expansion of k(u>) and u 2 
around cjq. By inverse Fourier transformation, we obtain the following wave equation in the 
time domain 





k^M 

n\ 



A(r, t) + 


1 

e 0 c 2 


where 


Defining the terms 


fc (n) M = 


dVcjtjj) 

doj n 


k {0) (u o) = 4 + iy, 
fc^( w o) = ki+i~, 


we obtain the operator 


D(u 0 ) = i 


.OiQ 


Oi 

2 


dt + J2 

n —2 


n! 


Using Eq. (17-19), the wave equation can be expressed as 



2 

B{?,t) = 0, 


(15) 


(16) 


(17) 

(18) 

(19) 


V 2 X + f — H- 2z/c 0 — k\-zr+iD 


'd ,9 ,A\ 
ai +kl m~' D i 


4 (f ' 4 )+ ^f'l +w “) B{r ' t)=0 - < 20 > 


[dz - - l dt 

In the moving reference frame defined by ,T = t-k x z, £ = z, the derivatives become 


_d ^ _ ,_d_ d__d_ 

dz ~ l dT' dt~dT' ( 21 ) 
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Then the wave equation becomes 


a 


'$]_+[ —+2ik 0 -2ki —+iD\ \--iD 


a 

dT 


d 


B(r ' T)=0 ' 


Reordering the terms 


... d ... 2k 0 d 2ko d n . d 

k ° - 2k 'df = " 2 


i 9 


a 


-2i,— _2i* 0 1 1 + — — 1+2 I --k. 


( 22 ) 


(23) 


and using Eq (23) in Eq (22), we obtain 


+ 2ik(\ [ 1 + 


i d 


a 


0 >odT \d£ 


-iD 


A(f,T) + 


OJn 


dt 


■f" iD -J- 2 


% 
, wo 


-h 


)A 1 

\l-iD 1 

) dT 

k J 


eoc 
A(r,t). 


1 + zw) B(f ' T) 


From the definition of k(u), we obtain 

dk 
duo 


w dyA yji 
c doj c 


where 


ki = Re 


k - ko 

K \ — —i 

W 0 


dsfi 


du 


UJ o' 


wo ko_ 

C UJq 


if 

OU CO 


If the right-hand side of Eq. (24) is very small, then Eq. (24) becomes 

-l 


<9,4 


iD + 


2k(] 


1 + 


_lJL 

w 0 aT. 


vi 


A + ±^L( 1 + ±i: 

2k 0 e 0 c 2 l u) 0 dT. 


B. 


If we assume 


(i + i±v 

V uiodTj 
A -tto , k 2 

D ~ t— + — 
2 2 


i-±± 

w 0 ar’ 

li 2 4 . h (iJL' 


,‘drj + 6 


dT, 


then Eq. (28) becomes 
dA 




a 0 M d 2 k 3 a 3 i 


2 dT 2 6 aT 3 2k, 


k 0 \ w 0 ary 




i U). 


i d 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 


75 1+-^? )B- (31) 


2A; 0 €q£ 2 \ cjq dT 
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3 NONLINEAR POLARIZATION 


If we consider only the third-order nonlinear effects governed by x (3 \ then the nonlinear 
polarization becomes 

Pni(f,t) = e of [ [ X {3) (t-ti,t-t 2 ,t-t 3 ):E{r,t 1 )E(r,t 2 )E(r,t z )dt 1 ,dt 2 ,dt 3 . (32) 

J—oo J —oo J -oo 

We further assume the following functional form for the third-order susceptibility, 

X i3) (t-h,t-t 2 ,t-t 3 ) =x {3) R{t-h)6(t-t 2 )5(t-t 3 ), (33) 

where 5(f) is the delta function and R(t) is the nonlinear response function normalized as 
J^oo R(t)dt = 1. Then the nonlinear polarization is given by 

Pni{r, t ) = 6oX (3) E(r, t) f R(t - t')E 2 (r, t')dt\ (34) 

J — OO 

here we assumed E || P n i- If the response function R(t) is nearly instantaneous, the we can 
apply the following Taylor series expansion 

R(t — ti)E 2 (r,ti)dti = f R(t — ti) 1 + dtiE 2 (r,t) 

J -°° [ n^i n ■ 9t n J 

~ ( 35 ) 

where T R = / 0 °° tR(t)dt. Then the nonlinear polarization becomes 

p ni{r, t ) « e 0 x {3) E(r, t) ^1 - T R ^J E 2 (f, t). (36) 

From the slowly varying envelope approximation, we obtain the following equations 

E 2 = ^4 2 e 2l ( fc oz-«;o t) _|_ 2 |j 4| 2 -I- (j4 *^ e -2i(koz-wot)^ (37) 

E^I-tJ^E 2 = |j_ ("1 _2 Mj ) j A 2 

+ A-{l-T K g-2io, 0 j J A 2 + 221 (f - T R ^j | A\ 2 

+ A jl - T r (^ + 2io,„j | {A*) 2 + 221* (l - T r P\ |2l| 2 

+ e- 3i (‘«*-“»‘)2l* |l - QL + 2io, 0 j j (A’) 2 . (39) 
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Assume that each c **(*o*-«o*) ( n = 3, l, - l, - 3 ) term can be separated from each other, 
then 


B ~ e 0 x^A 


r)A r) A* 

(3 + 2ioj 0 T R )\A\ 2 - 4 T r A* — - 2T r A-~- 


(40) 


If = 0, then B = 3e 0 X^M| 2j 4 considering to the Kerr nonlinearity. Using Eq. (40) in 
Eq. (20) we can obtain wave equation which expressed in terms of A(r,T ) only. 


4 SCALED NONLINEAR WAVE EQUATION 

For cylindrically symmetric system, we define the following dimensionless parameters 

£ r± 


V=T~, P 

J-'df To 


T n A 
r = —, Q = —, 

J-o -^o 


(41) 


where L df , r 0 , T 0 , A 0 are given by the input pulse and r ± = yfx 1 + y 2 . Using Eq. (41) the 
wave equation Eq. (31) can be rewritten as 


, AA.+_ 

df dr) ~\ 2 2T 0 2 dr 2 6T 0 3 dr 3 2k 0 \ u 0 T 0 dr) r 2 * p 


i d\ 1 _ 2 )^ , i (UqX ( 3) Mo| 


-^V 2 AQ4 


2/c 0 


B\ (42) 


and B' is given as 


i d 


B' « (l + AAlQ ( 3 + 2za;or ii )|Q| 2 -4^Q^-2^Q 


3Q 


to 0 To <9r y 
C 1 + -i(u 0 T/j)|(3| 2 + (^2^ 


AA.AnA 

To^ T 0 Q dr 




(43) 


where d 2 jdr 2 and higher order terms are neglected. Now we set L df = korl/2 and introduce 
the following parameters 


r TxoLdj , 'y r2 L d/ , $ tj oj T ' 


m 5 P n 7 0 (44) 


To’ 


fc 0 c 2 


Then Eq. (42) becomes 


d_Q 

dr] 


T iy d 2 5 d z if . d \ _ 

'2 “ ~2d? + 6d? + 4 ( 1-W7 ar) 


+ip {0 + lr) ie|2 +(tff - Tr) [ 2Q ' ^ + }«' < 45 > 
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Next we define two linear operators £>^(r) and Ddfij) and one nonlinear operator N(t) 
such as 


y. , x i'y d 2 8 d 3 T 

ds ^ ~ 2dr 2 + %d^~r 

(46) 

|) v ?. 

(47) 

N(t) « ip\Q\ 2 . 

(48) 

In the frequency domain, the linear operators become 


Afe(w) = —u + —u — —, 

(49) 

D d f(uj) = I(l-orw)Vj, 

(50) 


here we use F — (— iu>) n and u is a dimensionless angular frequency which correspond¬ 
ing the dimensionless time r. 


5 RSA RATE EQUATION 

Rate equations for a five-level RSA are 


dN 0 

dt 

— z No + kioNi + 

ricuo 

(51) 

dN x 

dt 

~ hu o N ° ( nu 0 + kl ° + kn)Nl + k2lN 2 ’ 

(52) 

dN 2 

dt 

= - k 2l N 2 , 

rtcjo 

(53) 

dN 3 

dt 

= ~ ^ Own ^' 30 ^ 3 k^N^, 

(54) 

dN 4 

dt 

= -Ns-k.M, 
hu)o 

(55) 


where Nj is the electron number density of the state j, a jk is the absorption cross-section 
for electron pumping from the state j to the state k and k jk is the decay rate from the state 
j to the state k. The assignment of the electron densities, Nj, in Fig. 1 is as follows: N 0 
corresponds to the ground level of the singlet state, S 0 ; N x corresponds to the first excited 
level of the singlet state, S \; iV 2 corresponds to the second excited level of the singlet state, 
S 2 ; N 3 corresponds to the first excited level of the triplet state, 7\; and JV 4 corresponds to 
the second excited level of the triplet state, J 2 . In these materials, the total number, N T , 




of electrons is conserved such that N 0 + M + JV 2 + N 3 + N 4 = N T . The initial population 
at 2 = 0 (i.e. before the pulse is incident on the RSA material) is given by N 0 = N T and 
Ni = N 2 = N 3 = N 4 = 0. The optical pulse has an input energy flux density I with angular 
frequency u 0 , h is Planck’s constant and I/huj Q is the input photon flux. 

The propagation equation of the energy flux density in the moving frame becomes 

dl 

= -(o'oi Nq + <7 12 iVi + g za N 4 )I. ( 56 ) 

If we assume / = ahu 0 \A\ 2 and a is a constant, then 


dA 1 

— ~ 2 ^ 01^0 + &12 Ni + g 34 N 4 )A. 

We can rewrite the rate equations and the propagation equation such as 

I 


MN = G + ~H N = [(? + a\A\ 2 fi] N, 


ON 
~dt 

dA 1, - 

ae = - 5 v-m, 


where 


N = 


and 


(No) 


0 

koi 

0 

kzo 

0 ) 


( ~coi 

0 

0 

0 

0\ 

Ni 


0 

— (kw + kiz) 

^21 

0 

0 


C’oi 

- 0\ 2 

0 

0 

0 

n 2 

> G — 

0 

0 

—k 2 i 

' 0 

0 

,H = 

0 

°\2 

0 

0 

0 

Nz 


0 

^13 

0 

—kzo 

&43 


0 

0 

0 

-<7 34 

0 

\ N *J 


[0 

0 

0 

0 

—k 43 ) 


^0 

0 

0 

°34 

oy 


& ' N — <Jq\Nq + (T\ 2 N\ + CT34./V3. 

Let’s introduce a following dimensionless parameter 

«. = #■ 

a t 


(57) 


(58) 

(59) 


( 60 ) 


(61) 


(62) 


Using dimensionless parameters defined in Eq. (41) and Eq. (62) the rate equations and the 
propagation equation become 


dr 

dQ _ L df N T 
drj 2 


= T 0 [G + cx\Ao\ 2 \Q\ 2 h} K, 
(d ■ K) Q. 


(63) 

(64) 
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Eq. (57) describes the propagation of a plane wave. However, for a three-dimensional pulse 
propagation in RSA material imbedded in a nonlinear medium, the pulse propagation equa¬ 
tion becomes 


dQ - , 

— = D ds + D d f + N + N rsa ] Q, 


(65) 


where N rsa — -a- #L d fN T /2 and the pulse propagation equation is coupled to the following 
rate equations 


dr 


= T 0 [g + c ,|4„| 2 |Q| 2 #]r 


( 66 ) 


6 SPLIT-STEP METHOD 


The nonlinear equation Eq. (65) can be expressed as 


dQ 

drj 


Dds + Ddf + N, 


'nl 


Q , 


where N ni = N + N rsa . The solution is described by 


(67) 


a Tj+Ar] A \ 

(Dds + D df + N nl )dT]') Q(ri) 

(Dds + D df ) At/ + j N nl dr) j Q(t/). (68) 

Generally linear operators D ds , D d j and nonlinear operator N n i do not commute. In order 
to solve the propagation equation, we use the symmetric split-step method given by 

Q(r) + A77) ^ e ^ Ddf e^~ DdsNnldr} ( 69 ) 

Linear operators D d f and D ds are solved in the frequency domain and the nonlinear operator 
N n i is solved in the time domain. The nonlinear operator N n i is solved using the trapezoid 
rule and then iterated until the solution converges 

r , , , An . „ 

Nni [Q(v )J dr, « Y \ N -i IQ(V)] + N nl [Q(rj + At/)]} . ( 70 ) 

In the first iteration N n i [Q(r] + At/)] is set equal to N n i [Q(t/)] to obtain an initial value for 
Q(r) + At/) which is then used for N n i [Q(t] + At/)]. The process is repeated until convergence 
is achieved. 
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7 NUMERICAL CALCULATION 

We consider a discrete time and space grid such as 


T->Ti, T]~^T] n , p -> pj. ( 71 ) 

For the fast Fourier transformation (FFT) algorithm, there must be exactly an integral power 
of 2 points equally spaced in the time domain 

Tj — iAi + T m i n . ( 72 ) 

The r] n and pj grids may be dynamically adjusted. 


7.1 Diffraction Operator 

The linear diffraction term in the the frequency domain is given by 

^ = xMVjQ(j7,p,w), 


( 73 ) 


where Q(r],p,u) is the Fourier transform of Q{t],p,t), x(w) = *(1 - <ju)/4 and Vj is the 
radial Laplacian with azimuthal symmetry 


V 2 = + — 

p p dp dp 2 ' 


( 74 ) 


This equation is solved using the following finite difference scheme. 

dQ _ QL +1 ~ % 

drj Ap 

and the Crank-Nicholson method 

V 2 p Q(r],p,u) = ^(V 2 p Q{r] n ,p,u) + V 2 p Q(r] n + Ar],p,co)), 

1 1 


( 75 ) 


( 76 ) 


V 2 pQ(Vn, P, u) - ^ 4Ap (Qlj +1 Qij-i) + 2{Ap) 2 (^ +1 “ + ^i-i) > ( 77 ) 


Pj 4 Ap 

where Ap = p J+1 — pj. Then 


_ Qij X(^i) r/vi+l An+1 . An An l 

Ar/ ” 4 p,Ap l ^' +1 

1 + Q?£, + < 3 ” i+ , - 24 "i + <&-,]. ( 78 ) 
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This equation can be rewritten as 

Uu-lQ&i + = V L-iQij-i + V iM + ^ + .< 3 ”j + .. ( 79 ) 

where 

xMAt? /_1_1_\ 

2Ap \2p, ApJ’ 

X(wi)Ap 


r/ 1 '. — 


u*. 

3,3 


rfi 

U 3,3+1 


VP , = 

JJ-1 


= 1 + 


v?. 

3,3 


= 1 - 


(Ap) 2 ’ 

Xfo)A t? /_1_1_\ 

2Ap V2 Pi ApJ 

x(^)At7 /_1_1_\ 

2Ap \2py Apy ’ 

x(wi)Ar/ 


V?. , 

V 3,3+l 


(80) 

(81) 

(82) 

(83) 

(84) 

(85) 


(Ap ) 2 ’ 

xfa) (J_ J_ N 

2Ap \2pj A p / 

which is valid for the interior points 0 < p { < p max . For on axis points (p = 0), the following 
facts are used; (1) for cylinderically symmetric systems dQ/dp = 0 at p = 0, and (2) the 
Maclaurin expansion for dQ/dp is used to define V 2 at p = 0 such as 

1 


1 


Q'(p) = Q , (o) + pQ"(o) + ^p 2 q"'(o) + 


lto-e'M = Q"( 0 ). 

Therefore, at the origin (p = 0), the cylindrical heat equation becomes 

^ = 2x{u) W 

Using the finite difference scheme described previously leads to 

ggzgu 2 xM 


Atj 


Pj+i—Pj-i 


Q’Sh-QZt'+QlM-Ql, ( QZJh-QS'+Qt^-Q. 

Pj +1 ~Pj Pj + i~pj 


i,3 


( 86 ) 

(87) 

( 88 ) 

• (89) 


At 0 = 0, i.e., j = 0 we can imagine extending the boundary to the left by a distance P-i 
which introduces the fictitious quantity Q"_,. This quantity can be eliminated using the 
fact that 

9 Q =0j <% - Qh- 1 , o, ( 3 ;, = (so) 

p -0 " l r-l 


dp 
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( 91 ) 


Therefore, if we assume that pi = p_i then the equation to be solved is 

Qto 1 - % _ 2 x(w 0 , 


Ar/ p\ 


, -[<3m 1 - 0“o 1 + 0"i - 0"ol- 
■'/ Pi 

Collecting terms this gives 

KoQlZ' + KM' = <0 %o + vio;,, 

where 

rri _ , . 2 X(^)Ap 

^0,0 — 1 + 9 7 

Pi 

TT i _ 2 X(^i)Ap 

^°,1 ~ p 2 > 

• T/i _ , 2 x(wi)A7] 

V 0,0 - 1 72 -> 

Pi 

2 x(^)At? 


T/ 2 — 

K 0,2 — 


Pi 


(92) 

(93) 

(94) 

(95) 

(96) 


For the maximum radial point (p = p^) the finite difference expression is derived using 
following procedures. To compensate for the extended boundary point p^+\ the following 
assumption are made. First, Pn+i is symmetrically placed around p N , i.e. 


Pn+i ~ Pn — Pn ~ Pn -ij Pn +i — 2pat — Pn-i- 


(97) 


Second, the value of the function Q" N - +1 is determined by linear interpolation through the 
points pM- 1 ; Pn and Pn+i, 


Q^n+i ~ Qj^n _ Qi^N ~ Qi,N -1 
Piv+i — Pn Pn ~ Pn -i 

Then the difference equation simplifies to 


Qi,N +1 — 2< 5r,/V “ Qi,N-\■ 


v'H.t'-iG&'-i+ u'k.nQ?? = 


^Tl + l _ 


where 





V i 

V N,N -1 




2 x(^)Ap 

2 Pn[Pn — Pn-\) ’ 

! _ 2 x(^»)Ar? 

2 Pn{Pn ~ Pn- i) ’ 
2x(^)Ap 
2 Pn{pn ~ Pn-i) ’ 

1 + 2x(^)Ap 

2Pn(Pn — Pn-i) 


(98) 

(99) 

( 100 ) 

( 101 ) 

( 102 ) 

( 103 ) 
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Now using the three finite difference expressions for the boundary points and the interior 
points, we can rewrite the expressions as a tri-diagonal matrix equation. 


m: 1 =v'qi, 


qt 1 = (uy^Qi 


(104) 


where matrices U and V 1 are tri-diagonal matrices which are determined by Eqs. (80-85), 
Eqs. (93-96), and Eqs. (100-103). Q " is a column vector whose elements are (• - • O”. . Q n . 
Qi,j +* / * Using a general matrix equation solver, we can easily solve the above equation. 


7.2 Dispersion Operator 

In the frequency domain, the dispersion operator is given by Eq. (50) and the differential 
equation for dispersion is given as 

dQ 

— = D ds (uj)Q. (105) 

The solution becomes 


Q(v + Ati/2, p, u>) = e^ u + s w3 2 ^Q(r],p,u>). (106) 

7.3 Nonlinear Operator 

The nonlinear operator is solved in the time domain using the trapezoid rule and then 
iterated until the algorithm converges. 

P+Ati „ „ 

J v dr) N\Q{rf)\ = -J-{N[Q( V )] + N[Q( V + Aiy)]}. (107) 

In the first iteration N[Q(ri+ A 77 )] is set equal to A^[<5( 7 ?)] to obtain an estimate for Q(r)+A t]) 
and the process is repeated until it converges. We consider a discrete time and* space grid, 

Q(r], p, t) -+ Q{ Vn , Pj , Ti ) = Q". (108) 

For N(r]), we write 

N(Vn) =ip\Q? j \ 2 . (109) 

The nonlinear part is comprised of simply multiplying the field in the time domain by the 
exponential 

JT V C(rj , r) = (w)^ r ) (n 0 ) 
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7.4 Rate Equation 

The rate equations are 


^ = T 0 {G + a\A 0 \ 2 \Q\ 7 H}K (111) 

To solve this equation, we consider the following grid in the r) — r plane where □ is the 
position of Q(77 n ,T t ) and O is the position of K(^ n+ i, ). Before the pulse enters the RSA 

material, all carriers are in the ground state. Therefore the initial condition for the carrier 
density is 




1 \ 
0 

0 

0 


( 112 ) 


VO ) 


For the propagation from Q(r} n ,Ti ) to Q(rj n ,T i+ 1 ), we assume the average carrier density is 
|[N(*7„-i, r i+ i)+K(j 7 n+ i, r i+ i)]. For the carrier excitation from N(7? n+ i, r f _i) to K(r) n+ i, r i+ i_), 

we assume the average field is \[Q{r] n , r i )+Q(r? n +i, t*)]- The values Q(^ n +i, a) and T i+\) 

are not determined before calculation. For the first iteration, <5(7?n+i, r ;) is set equal to 
Q(r/ n ,Ti). The rate equations and the pulse propagation equation are iterated until the 
solutions converge. From the rate equation, we can get 


H(7y n+ i,r i+ i) = exp \ To J* i+ ^ {G + a\A 0 \ 2 \Q\ 2 H} dr j N(7j n+ i,r i _i) 

= exp ^GToAr + Ha\A 0 \% |Q| 2 dr'j K(t 7n+ i,r i _i) 

« exp(GT 0 AT+ffa|/lo| !! rc,ATi{|Q(i)„,r i )| 2 +|(J(?)„ +I ,T i )P}) 

(113) 


00 I „ 


fc =0 


where R = GT 0 Ar + #c*|Ao| 2 ToAt| {|Q(?7„, Vi)| 2 + |<3(?7 „ + i, Tj)| 2 }. We should choose Ar 
such that 

|£|<1. (114) 

In our calculation, we keep terms up to k = 2. 
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8 ITERATION 

The rate equations and the wave equation are given by 


— To {G 4- a|j4o| 2 |Q| 2 tT} N, 

= [ D ds + D d f + N + N rsa ] Q , 


dQ 

dr] 


(115) 

(116) 


If we consider only the N rsa term in the wave equation, then 


9Q_ 1 c r n _ L df N T _ - 

Q r — NrsaQ — — cr • N Q. 


(117) 


The solutions are 

« exp(ToATG+ToArL 4 f ^{|Q(77 n ,t i )| 2 +|g(77 n+ 1 ,^)| 2 }^-%^,T^i), (118) 
Q(Vn+i,ti) « exp (-^a. {% n+ i, ri _i) + % n+ i,r i+ i)}) Q(r) n ,ti). (119) 

The k-th and (k+l)-th iterated solutions are 

t^j) »exp(T 0 ATG+r„ArI tf t{|Q(,„, (,)|N -Nfa*!,r^), (120) 

C <W) (>)»H,(i)«exp(-^.{% +J ,T j _i) + it(‘)( 17 „ + ,,r j+ i)})(3( % ,i i ). (121) 

Now consider the following errors 

£fy?„ + i,T i+ i) - K (fc) (j 7 n+ i,r i+ i) » 

T 0 ArL df ^ {|Q(JM+i,*i)| 2 ~ |<5 W (77 n+ 1 ,t i )| 2 } • N (fc) (7? n+ i, r^i), (122) 

IQ(?7n+i,ii)| 2 - |Q ( * +1) (»7n+i,ti)| 2 » 

—A 77 ^ ^ • IT - K(77 n+ i,Tj_i)|(5(r7 n ,ti)| ||Q(? 7 n+1 , tj )| 2 — |<2^(i7n+ii ^)| 2 } • (123) 


Then 


A _ _A?? 
AQ(*> ~ At 


2 a • • K(?7 n+ i, t^i)|Q(t7„, ^)| 2 , 


(124) 


where A<2 ( * +1 ) = |<5(77n+i,ti)| 2 -|Q (fc) (%+i,^)| 2 . If the right hand-side of Eq. (127) is smaller 
than 1, then the iteration improves the accuracy of the solution . 
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9 SOLITON SOLUTION 


To check the accuracy of our numerical method, we used the nonlinear Schrodinger equation 
(NLS) which has analytic solutions. By comparing the numerical solution to the analytic 
solution we verify the accuracy of our numerical method. The NLS equation is easily derived 
from our nonlinear wave equation and given by 


.dU 


i 


drj 


1 d 2 U 
+ 2 dr 2 


+ \U\ 2 U = 0. 


(125) 


For certain case, Eq. (128) solitons. The fundamental soliton is given as 


f 


U(t),t) — sech(T)exp(irj/2) 


(126) 


The intensity is independent of propagation distance rj. Figure 3 shows the analytic solution 
and the numerical solution of Eq. (128). In the numerical calculation, we used a time step 
of St = 0.05 and a propagation step of Srj = 0.000314 for a total propagation distance of 
3.14. The numerical solution agrees with the analytic solution yielding an error less than 
0.5% for -2 < r < 2 as shown in Figure 4. The error can be reduced by choosing a smaller 
propagation step and/or a smaller time step. In this calculation, we found that our numerical 
methods describe the soliton propagation very accurately. 


10 CONCLUSION 

In this report we described a numerical method for short optical pulse propagation in nonlin¬ 
ear media with RSA materials. Numerical solutions for soliton propagation are in excellent 
agreement with analytical solutions. In our numerical method, higher order group velocity 
dispersion terms, higher order nonlinear terms and the Raman term are easily incorporated 
for ultra-short pulse propagation. Five-level model reverse saturable absorbers(RSA) show 
optical limiting effect for appropriate input power intensities. For high power input pulse, 
the output pulse power is decreased significantly. Using our numerical results, we can design 
effective optical limiters for a given range of pulse powers. 
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Figure 1: Schematic energy level for a chromophore. The electronic states are represented by 
solid horizontal lines and the vibronic states are represented by dotted horizontal lines. Si 
represents a singlet state and Ti represents a triplet state. The photon absorption excitations 
are represented by solid vertical lines, and the decay processes are represented by wavy lines. 
The absorption coefficients cr^ and the decay constants kij are described in the text. The 
physical values used for our calculations are: cr 0 1 = 2.4 x 1()~ 18 cm 2 , a 12 = 3.0 x 10“ 17 cm 2 , 
<r 34 = 4.8 x 10~ 17 cm 2 , k 10 = 0.144ns" 1 , k 2l = l.Ops" 1 , k 13 = 77.8ms" 1 , k 30 = 50.0ms"\ 
^43 = 1-Ops \ and kn, k 22 , k 33 are due to vibrational decays which are assumed to be 
instantaneous. 
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Figure 2: Space-Time grid for pulse propagation. Squares represent coordinates where pulse 
is calculated and circles represent coordinates where carrier densities are calculated. 



Figure 3: Analytic and numerical solution for the fundamental soliton propagation. The 
solid line shows the analytic solution and the circles show the numerical solution. 
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